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^T" 1 Introduction 

^—^ In this article we study adaptive finite element methods (AFEM) with inexact solvers 

fH for a class of semilinear elliptic interface problems. We are particularly interested in 

^ nonlinear problems with discontinuous diffusion coefficients, such as the nonlinear 

S Poisson-Boltzmann equation and its regularizations. The algorithm we study con- 

, , sists of the standard SOLVE-ESTIMATE-M ARK-REFINE procedure common to 

many adaptive finite element algorithms, but where the SOLVE step involves only a 
^-H full solve on the coarsest level, and the remaining levels involve only single Newton 

^ updates to the previous approximate solution. We summarize a recently developed 

AFEM convergence theory for inexact solvers appearing in [2], and present a se- 
quence of numerical experiments that give evidence that the theory does in fact pre- 
fvq diet the contraction properties of AFEM with inexact solvers. The various routines 

1^ used are all designed to maintain a linear-time computational complexity. 

^^ An outline of the paper is as follows. In Section 2, we give a brief overview of 

^-H the Poisson-Boltzmann equation. In Section 3, we describe AFEM algorithms, and 

^■H introduce a variation involving inexact solvers. In Section 4, we give a sequence of 

^ numerical experiments that support the theoretical statements on convergence and 

optimality. Finally, in Section 5 we make some final observations. 



CO 



t 



2 Regularized Poisson-Boltzmann Equation 

We use standard notation for Sobolev spaces. In particular, we denote || • ||o,g the L^ 
norm on any subset G C M^, and denote || • || i,2,g the H^ norm on G. 

Let Q := Q^ UP U Qs be a bounded Lipschitz domain in M?, which consists 
of the molecular region 12^, the solvent region Qs and their interface F := Qm^^s 
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Fig. 1. Schematic of a molecular domain. 



(see Figure 2). Our interest in this paper is to solve the following regularized Poisson- 
Boltzmann equation in the weak form: find u e H^{Q) := {u e H^{Q) : u\^q = g} 
such that 

^(w,v) + (%),v) = (/,v) \/veH^{Q), (1) 

where a{u^v) = Jq eVw • VvJx, {b{u)^v) = Jq k^ smh{u)vdx. Here we assume that 
the diffusion coefficient £ is piecewise positive constant e|x2,„ = e^ and e|x2^ = e^. The 
modified Debye-Huckel parameter K^ is also piecewise constant with k'^{x)\q^ = 
and k'^{x)\q^ > 0. The equation (1) arises from several regularization schemes (cf. 
[4, 5]) of the nonlinear Poisson-Boltzmann equation: 



-V • (eVu) + K^ sinhw = ^Zi5{xi), 

i=l 

where the right hand side represents N fixed points with charges Zi at positions xt, 
and 5 is the Dirac delta distribution. 

It is easy to verify that the bilinear form in (1) satisfies: 

<^o||w||i,2 <a(u,u), a(u,v) < ci||w||i,2||v||i,2, ^u,veHl{Q.), 

where < cq < ci < oo are constants depending only on e. These properties imply 
the norm on //J (12) is equivalent to the energy norm ||| • ||| : //q (12) ^ R, 

infill =a(w,w), ^oll^lll 2 — infill ^<^l||^||l2- 

Let ^ be a shape-regular conforming triangulation of ^, and let yg(^) := {v G 
if ^ (^) : v|^ G Pi (t) Vt G ^} be the standard piecewise linear finite element space 
defined on 3h. For simplicity, we assume that the interface F is resolved by 3h. Then 
the finite element approximation of (1) reads: find Uh^yg{3h) such that 

a{uh,v)^(b{uH),v) = {f,v), VvGyo(5J,). (2) 

We close this section with a summary of a priori LT bounds for the solution u 
to (1) and the discrete solution Uh to (2), which play a key role in the finite element 
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error analysis of (2) and adaptive algorithms. For interested reader, we refer to [4, 7] 
for details. 

Theorem 1. There exist w+,w_ G U°{Q) such that the solution u of (1) satisfies the 
following a priori LT bounds: 

U-<u<u^^ a.e.inQ. (3) 

Moreover, if the triangulation ^ satisfies that 

a{(i>u(i>j)<--^ Y, 1^1, for some a > 0, (4) 

for all the adjacent vertices i ^ j with the basis function (\>i and 0^, then the discrete 
solution Uh of (2) also has the a priori LT bound 

hh\\L-(Q) <C, (5) 

where C is a constant independent ofh. 

We note that the mesh condition is generally not needed practically, and in fact can 
also be avoided in analysis for certain nonlinearites [2] . 



3 Adaptive FEM with Inexact Solvers 

Given a discrete solution w/^ G yg(^), let us define the residual based error indicator 

ri{uh,T): 

r]\uh,T)=h',\\b{uh)-f\\l,^ X he\\[{£^uh)-ne]\\l,, 

ecdT 

where [{sVuh) • ne] denote the jump of the flux across a face e of T. For any subset 
^ C 5/j, we set if'^Uh^y) := Y^TeS^V^i^h^^)- By using the a priori L°° bounds 
Theorem 1, we can show (cf. [7]) that the error indicator satisfies: 

|||w-w/,|p <Ciri^{uh,^h)'. (6) 

and 

|r7(v,T)-r7(w,T)|<C2|||v-w|||co„ Wv.weVgi^h) (7) 

where co^ = U^'e^^^f'nfy^^d"^' and |||v|||^^ = /^^ e|VvpJx. 

Given an initial triangulation ^, the standard adaptive finite element method 
(AFEM) generates a sequence [uk^ ^, {t] (wy^, t)}^^^^] based on the iteration of the 
form: 

SOLVE ^ ESTIMATE ^ MARK ^ REFINE. 

Here the SOLVE subroutine is usually assumed to be exact, namely u^ is the exact 
solution to the nonlinear equation (2); the ESTIMATE routine computes the element- 
wise residual indicator r7(wy^,T); the MARK routine uses standard Dorfler marking 
(cf. [6]) where ^k C ^ is chosen so that 



4 Michael Hoist, Ryan Szypowski, and Yunrong Zhu 

n(uk,^k) >Ori{uk,^k) 

for some parameter 6 G (0,1]; finally, the routine REFINE subdivide the marked 
elements and possibly some neighboring elements in certain way such that the new 
triangulation preserves shape-regularity and conformity. 

During last decade, a lot of theoretical work has been done to show the conver- 
gence of the AFEM with exact solver (see [9] and the references cited therein for 
linear PDE case, and [8] for nonlinear PDE case). To the best of the authors knowl- 
edge, there are only a couple of convergence results of AFEM for symmetric linear 
elliptic equations (cf. [10, 1]) which take the numerical error into account. To distinct 
with the exact solver case, we use Uk and ^ to denote the numerical approximation 
to (2) and the triangulation obtained from the adaptive refinement using the inexact 
solutions. 

Due to the page limitation, we only state the main convergence result of the 
AFEM with inexact solver for solving (1) below. More detailed analysis and exten- 
sion are reported in [2]. 

Theorem 2. Let {^,%}yt>o be the sequence of meshes and approximate solutions 
computed by the AFEM algorithm. Let u denote the exact solution and u^ denote 
the exact discrete solutions on the meshes ^. Then, there exist constants jl > 0, 
V G (0, 1), 7 > 0, and a G (0, 1) such that if the inexact solutions satisfy 

l^hk - uk\f + \\\uk^i - %+i IIP < vr7^(%, ^k) (8) 

then 

\\\u- uk^i f + r^^(%+i , ^m) < a^{\\\u- ukf + r^^(%, ^k))' (9) 

Consequently, limy^^oo u^ = limy^^oo Uk = u. 

The proof of this theorem is based on the upper bound (6) of the exact solution, 
the Lipschitz property (7) of the error indicator, Dorfler marking, and the following 
quasi-orthogonality between the exact solutions: 

|||W — Wy^+l||| < A|||t/ — Wy^ll — |||Wy^+l — Wy^ll (10) 

where A can be made close to 1 by refinement. For a proof of the inequality (10), see 
for example [7] . 

To achieve the optimal computational complexity, we should avoid solving the 
nonlinear system (2) as much as we could. The two-grid algorithm [11] shows that a 
nonlinear solver on a coarse grid combined with a Newton update on the fine grid still 
yield quasi-optimal approximation. Motivated by this idea, we propose the follow- 
ing AFEM algorithm with inexact solver, which contains only one nonlinear solver 
on the coarsest grid, and Newton updates on each follow-up steps: In Algorithm 1, 
the N SOLVE routine is used only on the coarsest mesh and is implemented using 
Newton's method run to certain convergence tolerance. For the rest of the solutions, 
a single step of Newton's method is used to update the previous approximation. That 
is, UPDATE computes %+i such that 



Inexact AFEM for Nonlinear PBE 



Algorithm 1 : %,^,{r7(%,T)}^^^ := lnexact_AFEM(^,0) 



1 uo = uo:= NSOLVE(^) %Nonlinear solver on initial triangulation 

2 for ^ :=0, !,••• do 

3 {r7(%,T)}^^^^ := ESTIMATE(%,^^) 

4 J^k := MARK({r7(%, t)}^^j. , #^, 0) 

5 #^+i:=REFINE(#^,^^) 

6 %+i := UPDATE(%, ^+i) %One-step Newton update 

7 end 



a{uk+\-Uk,(i>)^(b'(uk)(uk+\-Uk),(i>) =0 (11) 

for every G V (^+i ). We remark that since (1 1) is only a linear problem, we could 
use the local multilevel method to solve it in (near) optimal complexity (cf. [3]). 
Therefore, the overall computational complexity of the Algorithm 1 is nearly opti- 
mal. 

We should point out that it is not obvioius how to enforce the required approxima- 
tion property (8) that Uk must satisfy for the theorem. This is examined in more detail 
in [2] . However, numerical evidence in the following section shows Algorithm 1 is an 
efficient algorithm, and the results matches the ones from AFEM with exact solver. 



4 Numerical Experiments 

In this section we present some numerical experiments to illustrate the result in The- 
orem 2, implemented with FETK. The software utilizes the standard piecewise-linear 
finite element space for discretizing (1). Algorithm 1 is implemented with care taken 
to guarantee that each of the steps runs in linear time relative to the number of ver- 
tices in the mesh. The linear solver used is Conjugate Gradients preconditioned by 
diagonal scaling. The estimator is computed using a high-order quadrature rule, and, 
as mentioned above, the marking strategy is Dorfler marking where the estimated 
error have been binned to maintain linear complexity while still marking the ele- 
ments with the largest error. Finally, the refinement is longest edge bisection, with 
refinement outside of the marked set to maintain conformity of the mesh. 

We present three sets of results in order to explore the effects of the inexact 
solver in multiple contexts. For all problems, we present a convergence plot using 
both inexact and exact solvers (including a reference line of order N~ 3 ) as well as 
a representative cut-away of a mesh with around 30,000 vertices. The exact discrete 
solution is computed using the standard AFEM algorithm where the solution on each 
mesh is computed by allowing Newton's method to continue running to convergence 
to within a tolerance of 10~^. This modifies not only the solution on a given mesh, 
but also the sequence of meshes generated, since the algorithm may mark differ- 
ent simplices. However, in all cases, convergence is identical to high precision and 
meshes place the unknowns as expected. 
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The first result uses constant coefficients across the entire domain Q = [0, 1]^ 
and we choose a right hand side so that the derivative of the exact solution is large 
near the origin. The boundary conditions chosen for this problem are homoegenous 
Dirichlet boundary conditions. Specifically, the exact solution is given by u = u\U2 
where 

u\ = sin(;rx) sin (ttj) sin(;rz) 

is chosen to satisfy the boundary condition and 
The results can be seen in Figure 2. 








Fig. 2. Convergence plot and mesh cut-away for the corner singularity problem. 



The second result uses the domain Q. = [—1,1]^ and £2^ = [—\^\] with con- 
stants £s = ^0^£jn = 2^Ks = 1, and Kfn = 0. Homogcneous Neumann conditions are 
chosen for the boundary and the right hand side is simplified to a constant. Because 
an exact solution is unavailable for this (and the following) problem, the error is 
computed by comparing to a discrete solution on a mesh with around 10 times the 
number of vertices as the finest mesh used in the adaptive algorithm. Figure 3 shows 
the results for this problem. As can be seen the refinement favors the interface and 
the inexact and exact solvers perform as expected. 

The final result is chosen to test the robustness of the inexact method to large 
coefficient jumps. The domain and boundary conditions are the same as the previous 
example, but the constants chosen as are £s = 1000, £m = 10, k*^ = 1, and k^ = 0. The 
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Fig. 3. Convergence plot and mesh cut-away for the Poisson-Boltzmann problem. 



results can be seen in Figure 4, and they show a scenario very similar to the previous 
example, with the refinement restricted even further to the interface. 




Fig. 4. Convergence plot and mesh cut-away for the exaggerated-jump Poisson-Boltzmann 
problem. 
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5 Conclusion 

In this article we have studied AFEM with inexact solvers for a class of semilinear 
elliptic interface problems with discontinuous diffusion coefficients, with emphasis 
on the nonlinear Poisson-Boltzmann equation. The algorithm we studied consisted 
of the standard SOLVE-ESTIMATE-M ARK-REFINE procedure common to many 
adaptive finite element algorithms, but where the SOLVE step involves only a full 
solve on the coarsest level, and the remaining levels involve only single Newton up- 
dates to the previous approximate solution. The various routines used are all designed 
to maintain a linear-time computational complexity. Our numerical results indicate 
that the recently developed AFEM convergence theory for inexact solvers in [2] does 
predict the actual behavior of the methods. 
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